Helmholtz Problem

from netgen.occ import *
from ngsolve import *
from ngsolve.webgui import Draw
from ngsolve.bem import *

keys: Acoustic Scattering, Brakhage-Werner formulation, Helmholtz equation, Dirichlet data

Helmholtz Problem#

We consider a Helmholtz problem with a plane wave as initial condition, i.e., \(m = e^{i\,\kappa \,x} \) and

\[\begin{split} \left\{ \begin{array}{rcl l} \Delta u + \kappa^2 \, u &=&0, \quad &\Omega \subset \mathbb R^3\,,\\ \gamma_0\, u &=& m, \quad &\Gamma = \partial \Omega\, .\end{array} \right. \end{split}\]

Combined field integral equations combine single and double layer integral operators, one simple option is the Brakhage-Werner formulation.

The solution is represented as

\[ u = i \,\kappa \,\mathrm{SL}(j) - \mathrm{DL}(j)\,, \]

where \(j\) solve the boundary integral equation $\( \left( \frac{1}{2} + \mathrm K + i \, \kappa \,\mathrm V \right) \, \mathrm j = \mathrm u_{in} \qquad \text{on} \, \Gamma = \partial \Omega \)$

The single layer and double layer operator are defined using the fundamental solution of the Helmholtz operator,

\[ G(x,y) = \frac{e^{i\kappa \|x-y\|}}{4\pi\,\|x-y\|}. \]

The layer potentials \( \mathrm SL\) and \(\mathrm DL\) have the same form as in the Laplace case, except that the Laplace kernel is replaced by the Helmholtz kernel \(G\).

kappa=10
order=4
screen = WorkPlane(Axes( (0,0,0), Z, X)).RectangleC(15,15).Face()

sp = Fuse(Sphere( (0,0,0), pi).faces)
screen.faces.name="screen"
sp.faces.name="sphere"
shape = Compound([screen,sp])

mesh = shape.GenerateMesh(maxh=5/kappa).Curve(order)
Draw (mesh);
fes_sphere = Compress(SurfaceL2(mesh, 
                                order=order, 
                                complex=True, 
                                definedon=mesh.Boundaries("sphere")))
u,v = fes_sphere.TnT()
fes_screen = Compress(SurfaceL2(mesh, order=order, 
                                dual_mapping=True, 
                                complex=True, 
                                definedon=mesh.Boundaries("screen")))
print ("ndof_sphere = ", fes_sphere.ndof, "ndof_screen =", fes_screen.ndof)
ndof_sphere =  17010 ndof_screen = 31110
with TaskManager():
    C = HelmholtzCF(u*ds("sphere"), kappa)*v*ds
    u,v  = fes_sphere.TnT()
    Id = BilinearForm(u*v*ds).Assemble()
lhs = 0.5 * Id.mat + C.mat
m = exp(1j * kappa * x) 
rhs = LinearForm(m*v*ds).Assemble()
gfu = GridFunction(fes_sphere)
pre = BilinearForm(u*v*ds, diagonal=True).Assemble().mat.Inverse()
with TaskManager():
    gfu.vec[:] = solvers.GMRes(A=lhs, b=rhs.vec, pre=pre, maxsteps=40, tol=1e-8)
GMRes iteration 1, residual = 60.52440579389444     
GMRes iteration 2, residual = 24.00171041505929     
GMRes iteration 3, residual = 10.816186639702561     
GMRes iteration 4, residual = 5.178501301182413     
GMRes iteration 5, residual = 2.549268448878952     
GMRes iteration 6, residual = 1.2888054956882669     
GMRes iteration 7, residual = 0.6549832602491108     
GMRes iteration 8, residual = 0.3554164558330167     
GMRes iteration 9, residual = 0.19990866606606375     
GMRes iteration 10, residual = 0.10836869993844714     
GMRes iteration 11, residual = 0.05931445314039659     
GMRes iteration 12, residual = 0.03294608123183362     
GMRes iteration 13, residual = 0.01381913910824682     
GMRes iteration 14, residual = 0.007360642470274796     
GMRes iteration 15, residual = 0.0033138105281745656     
GMRes iteration 16, residual = 0.0019472886266694597     
GMRes iteration 17, residual = 0.0011529608192671487     
GMRes iteration 18, residual = 0.0006865451125338563     
GMRes iteration 19, residual = 0.00040647668374399725     
GMRes iteration 20, residual = 0.00023525302450811214     
GMRes iteration 21, residual = 0.0001256833512855774     
GMRes iteration 22, residual = 7.272844838705902e-05     
GMRes iteration 23, residual = 4.282868528421695e-05     
GMRes iteration 24, residual = 2.238947095910066e-05     
GMRes iteration 25, residual = 1.132365408107729e-05     
GMRes iteration 26, residual = 5.681772319703246e-06     
GMRes iteration 27, residual = 2.8154812159697885e-06     
GMRes iteration 28, residual = 1.3918984915180512e-06     
GMRes iteration 29, residual = 7.931632776571544e-07     
GMRes iteration 30, residual = 4.523061248529777e-07     
GMRes iteration 31, residual = 2.723734968721736e-07     
GMRes iteration 32, residual = 1.609611766851694e-07     
GMRes iteration 33, residual = 9.476709863742893e-08     
GMRes iteration 34, residual = 5.2145026155061434e-08     
GMRes iteration 35, residual = 2.970933842689033e-08     
GMRes iteration 36, residual = 1.7339940734545346e-08     
GMRes iteration 37, residual = 9.167410459416053e-09     
Draw (gfu, order=5, min=-1, max=1);

Postprocessing on screen#

uscat = GridFunction(fes_screen)
with TaskManager():
    uscat.Set(1j*kappa*HelmholtzSL(u*ds("sphere"),kappa)(gfu) - HelmholtzDL(u*ds("sphere"), kappa)(gfu), \
              definedon=mesh.Boundaries("screen"))
print ("Scattered field")
Draw (uscat, mesh, min=-1,max=1, animate_complex=True, order=4);
Scattered field
uin = mesh.BoundaryCF( {"screen": m }, default=0)
print ("Total field")
Draw (uin+uscat, mesh, min=-1,max=1, animate_complex=True, order=4);
Total field